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We present the first numerical solution to the next to leading order Balitsky-Kovchegov (BK) 
equation in coordinate space in the large-W limit. In addition to the dipole operator we also 
solve the evolution of the “conformal dipole” for which the conformal invariance breaking double 
logarithmic term is absent from the evolution equation. The NLO corrections are shown to slow down 
the evolution. We show that the solution depends strongly on the details of the initial condition, 
and that the solution to the equation is not positive definite with all initial conditions relevant for 
phenomenological applications. 
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I. INTRODUCTION 

At the high energies of present day collider experi¬ 
ments, the available phase space for gluon radiation is 
very large. Since each emitted gluon is itself a source of 
further emissions, a typical scattering event involves an 
exponentially growing cascade of gluons. At high enough 
energy this cascade can fill up the available phase space 
to the extent that the gluons begin to reinteract. This 
is the origin of gluon saturation, where the phase space 
below an energy dependent transverse momentum scale 
Qs is dominated by nonlinear gluon interactions. This 
recombination restores the unitarity of the scattering S- 
matrix, which would be violated by an unlimited expo¬ 
nential growth of the cascade. 

At high energy, it is convenient to describe QCD scat¬ 
tering off a hadronic target using the eikonal approxi¬ 
mation. The natural degrees of freedom describing the 
gluons in the cascade are then transverse coordinate de¬ 
pendent Wilson lines which describe the eikonal prop¬ 
agation of a high energy quark through the target color 
field. Cross sections can be expressed in terms of correla¬ 
tors of these Wilson lines. These correlators are universal 
objects that enable a consistent description of many dif¬ 
ferent scattering phenomena, such as deep inelastic scat¬ 
tering (DIS) [1], single inclusive particle production in 
proton-nucleus collisions [2-5] and two-particle correla¬ 
tions [6-9]. In the dilute limit, when nonlinearities are 
unimportant, these correlators reduce to unintegrated 
gluon distributions [10]. This framework is often referred 
to as the Color Glass Condensate (CGC, for a review, see 
e.g. Ref. [11]). 

The dependence of these Wilson line correlators on en¬ 
ergy, or equivalently Bjorken x or rapidity y = In 1/a;, can 
be calculated perturbatively even in the nonlinear satura¬ 
tion regime. This energy dependence is described by the 
Balitsky hierarchy of evolution equations, which for many 
practical applications can be replaced by its mean field 
version, the Balitsky-Kovchegov (BK) evolution equa¬ 
tion [12, 13]. It resums large logarithms ~ as In 1/a; to 
all orders. Current phenomenological works typically use 


the leading order BK equation with the running coupling 
corrections derived in Ref. [14]. 

While leading order calculations can give a good phys¬ 
ical description of the process, higher order corrections 
can be numerically very large. It is therefore extremely 
important for quantitative comparisons with data to per¬ 
form the CGC calculations at next-to-leading (NLO) or¬ 
der accuracy in the QCD coupling ag. First steps in this 
direction have been taken in particle production in pA 
collisions [15-18] and DIS [19, 20]. 

A crucial ingredient of a consistent NLO treatment in 
this context is solving the NLO BK equation, which de¬ 
scribes the energy or rapidity, dependence of the Wil¬ 
son line correlators. This equation has been derived in 
Ref. [21]. Its linearized limit, the NLO BFKL equation 
has been known already before [22-24]. It has been noted 
that it is affected with contributions from large transverse 
momentum logarithms, a feature that seems to remain 
true when the NLO BFKL equation is complemented 
with an absorptive boundary condition to emulate the 
nonlinear effects [25]. Elaborate resummation schemes 
have been proposed to treat these contributions [26-29]. 
These resummations, however, rely on the linear struc¬ 
ture of the equation and cannot easily be generalized to 
BK. There have also been proposals for a kinematically 
constrained version of the BK equation [30, 31] to solve 
these issues. While it is expected that these same in¬ 
stabilities also manifest themselves in the fully nonlinear 
BK equation in coordinate space, this has never been 
shown explicitly. In order to directly address this ques¬ 
tion, we will in this work solve numerically the next to 
leading order BK evolution equation in the form derived 
in Ref. [21]. We will also study numerically the equa¬ 
tion for the “composite conformal dipole”, proposed in 
Ref. [32] to address some of the issues with the NLO BK 
equation. The precise forms of the equations studied nu¬ 
merically are written down in Secs. II and III. Our nu¬ 
merical results are then discussed in Sec. IV. In princi¬ 
ple the dipole scattering amplitude obtained here could 
be convoluted with the NLO photon impact factor from 
[19, 20] for a comparison with experimental data. We 


2 


will demonstrate, however, that this would be difficult 
since the equation itself is unstable for many values of 
the initial conditions in the phenomenologically relevant 
regime. 


Here we have taken the large-A^c and mean field limits 
to express the equation in a closed form in terms of only 
the correlator of two Wilson lines (the “dipole”): 


II. BALITSKY-KOVCHEGOV EQUATION AT 
NLO 


We study in this work the next-to-leading order BK 
evolution equation derived in [ 21 ], which we write as: 


S{r) = ^^{TTU,Ul}. 


( 2 ) 


fy /V 

dyS{r) = ® [S{X)S[Y) - S{r)] 

a^N ^ 

+ 0 [^(X)5(z - z')S{Y') - S{X)S{Y)] 

+ “ 1 ^^/ ® - ^(^)] ( 1 ) 


Here the brackets () stand for an average over the target 
color field. The kernels appearing in Eq. (1) are 
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The coordinates are two dimensional vectors denoted as 
r = x — y,X=x — z,Y = y — z, X' = x — z' and Y' = y — 
z'. The convolutions ® are calculated by integrating over 
the vectors z and z'. The kernel Ki consists of the leading 
order BK kernel r^/(^ 2 ^ 2 ^ 3 ^^ NLO correction ~ as, 

and the beta function coefficient is /3 = yA"c — |?T-f with 
nf = 3. 

Part of the NLO corrections, especially the term in¬ 
volving the renormalization scale fj?, should be absorbed 


into the running of the strong coupling ag. What terms 
exactly are absorbed in Og, and at which scale it is evalu¬ 
ated, is a scheme choice. We adapt the choice derived in 
Ref. [14] and replace all terms in Ki proportional to the 
(3 function by the Balitsky running coupling prescription 
that is also used to solve the leading order BK equation 
with running coupling corrections. For the other terms 
we choose to evaluate Og at the scale given by the size of 
the parent dipole r. We thus write the first kernel as 
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We use the same expression of the coupling in terms write in terms of S, the results of our calculation will 
of as in Ref. [33]. While the equation is simpler to 
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Figure 1: Dipole amplitude and conformal dipole amplitude 
at initial condition and after evolution compared to the solu¬ 
tion of the LO BK equation at rapidities 1 / = 0, 5 and y = 10 
(from right to left). 


be expressed in terms of the scattering amplitude^ N = 
1-S. 

At finite Nc, correlators of up to six Wilson lines would 
also be needed in order to evaluate the rapidity derivative 
of the dipole operator in Eq. (1). In principle one could 
obtain the higher-point functions from a solution of the 
NLO JIMWLK equation [34, 35], or using e.g. a Gaussian 
approximation which allows one to write any higher point 
function in terms of the dipole operators as described in 
Ref. [36]. The contribution from finite-A), corrections 
to the leading order BK equation have been studied in 
Ref. [37], and their contribution is numerically found to 
be even smaller than the ~ 1/A^c^ one would naively 
expect, so we feel justified in neglecting them here. 


III. EVOLUTION EQUATION FOR THE 
CONFORMAL DIPOLE 

The Wilson lines are conformally invariant, and thus 
their evolution equation should be conformal in a con¬ 
formal field theory. In QCD, one expects the conformal 
invariance to be broken only by the running of the cou¬ 
pling. However, the evolution Eq. (1) also has a confor¬ 
mal symmetry breaking NLO double logarithmic term 
InX^/r^ InY^lr^ in the kernel Ki. Diagrammatically 
this contribution arises from the diagrams with a loop 
in the 1 —>■ 2 dipole transition where one gluon interacts 
with the shockwave, see discussion and Eig. 9 in Ref. [21]. 


^ Note that a term N(Y') — N{Y) is missing from Eq. (136) of 
Ref. [21] 


Q.!.o/Aqcd = 19,7 = 0.6 



Figure 2: Evolution speed for the conformal and non- 

conformal dipoles as a a function of the saturation scale com¬ 
pared to the leading order BK equation solution. 

The reason for the conformal invariance breaking is the 
fact that the derivation of Ref. [21] uses a cutoff in the 
longitudinal direction that violates the symmetry. This 
was confirmed by the appearence of the same nonconfor- 
mal double logarithm in the fully conformally invariant 
V = 4 supersymmetric Yang-Mills (SYM) theory [32]. 

A possible way to restore the conformal invariance, 
proposed in Ref. [32], is to rewrite the evolution equa¬ 
tion in terms of the conformal dipole defined as 

Here a is an arbitrary dimensional constant which will 
eventually cancel from the evolution equation. Using 
Eq. (1) one can then derive the NLO evolution equa¬ 
tion for the conformal dipole. The resulting equation 
turns out to differ from the NLO BK equation only 
by the disappearance of the double logarithmic term 
InX^/r^ In Y^/r^ from Ki, and the appearance of an ad¬ 
ditional contribution 

2r2 r^(z-z')2 

X2Y'2(z- X'2Y2 

in the kernel 7^2 • Now the only term that breaks the 
conformal invariance is the running coupling as- The 
corresponding evolution equation in A = 4 SYM theory 
is fully conformal [32]. 

IV. NUMERICAL ANALYSIS 

We solve the evolution equations for the non-conformal 
and conformal dipoles on a logarithmical grid in r using a 
















4 





Figure 3: Logarithmic derivative of the dipole amplitude (evolution speed) at initial condition with different values for the 
anomalous dimension. 





Figure 4: Evolution speed of the dipole amplitude at y = 5 with different values for the anomalous dimension at the initial 
condition. 


Runge-Kutta method. The four-dimensional integral in 
the NLO part is computed using an adaptive Monte Carlo 
algorithm. As an initial condition we use a McLerran- 
Venugopalan model [38] 



modified by introducing an anomalous dimension 7 which 
controls the power-like tail of the dipole amplitude for 
small dipoles. This parametrization is used in phe¬ 
nomenological hts to DIS data e.g. in Ref. [Ij. Deter¬ 
mining the correct values for q and 7 would require a 
full NLO fit to e.g. DIS data, which we are not perform¬ 
ing here. It is not obvious that the initial condition would 
be the same as for the leading order equation. We shall 
here merely explore the general behavior of the equation 
with different values for and 7 without aiming for 
phenomenologically relevant values in this work. 

We find that for some initial conditions (see discus¬ 
sion later) the evolution becomes unstable, such that the 
dipole amplitude starts to decrease and may even turn 
negative for small dipoles. It however follows from the 
definition of the dipole operator, Eq. (2), that one should 


have N(r) —>■0 in the limit r —>■ 0 , which is violated by 
non-zero amplitude at small r. Also the convolution with 
the kernel Ki in Eq. (1) does not converge if this require¬ 
ment is not fulfilled. To obtain this property, we freeze 
N{r) =0 in the region where the evolution would turn 
it negative. 

The dipole amplitudes from the NLO equation are 
compared with the solution of the leading order equa¬ 
tion in Fig. 1. The main effect of the evolution is to 
increase the amplitude at small r, while maintaining it 
below the black disk limit of = 1 at large r. This 
leads to the curve N{r) moving towards the left (smaller 
r) with rapidity in Fig. 1. It can be seen that the NLO 
corrections reduce the evolution speed significantly but 
the shape of the dipole amplitude remains roughly un¬ 
changed. The solution in Fig. 1 has an initial condition 
Qs.o/Aqcd ~ 19, 7 = 0.6, deliberately chosen such that 
the dipole amplitude increases at small dipoles through¬ 
out the evolution over the rapidity interval studied here. 
The evolution for the conformal dipole with the same ini¬ 
tial condition is also shown. Note that at leading order 
also the conformal dipole evolution is given by the stan¬ 
dard LO BK equation. In the LO BK equation, we also 
use the same “Balitsky” running coupling prescription. 
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Figure 5: Evolution speed of the conformal dipole amplitude at initial condition with different values for the anomalous 
dimension. 





Figure 6 : Evolution speed of the conformal dipole amplitude at y = 5 with different values for the anomalous dimension at the 
initial condition. 


The change of the saturation scale with energy is quanti¬ 
fied more precisely in Fig. 2 with the evolution speed of 
the saturation scale 


where the precise definition of used here is 

iV(r2 = 2/g2) = 1 _ e-i/2, (11) 

The NLO corrections can again be seen to significantly 
slow down the evolution speed. The conformal and “non- 
conformal” dipoles have comparable evolution speeds af¬ 
ter a few evolution steps, and the total evolution speed 
decreases slowly as a function of Qs- Note that the small 
anomalous dimension in the initial condition makes the 
leading order evolution faster than A ^ 0.2... 0.3 ob¬ 
tained in leading order fits with 7 ^ 1 [1, 3, 39, 40]. 
Also the parameter Q^q that controls the initial satu¬ 
ration scale is not the same as the saturation scale 
obtained by solving the equation ( 11 ), and in this case 
Qs.o/^qcd ~ 19 corresponds to having an initial satura¬ 
tion scale Qs/Aqcd 40. 

One would generally expect N to increase with rapid¬ 
ity, corresponding to the physical picture of more gluons 


being emitted when the available phase space increases 
with increasing collision energy. This is the behavior 
seen in the LO equation. To study when exactly this 
happens we show in Fig. 3 the evolution speed (logarith¬ 
mic derivative of the dipole amplitude dyN{r)/N{r)) at 
y = 0 with different values for the anomalous dimension 
7 and initial saturation scale Qs,o as a function of the 
parent dipole size. We see that the scattering amplitude 
does indeed increase, but only for a suitable choice of 
the initial conditions: small enough 7 and large enough 
Qsfl- Let us discuss the interpretation of the logarith¬ 
mic derivative plots in more detail. For smaller Qg the 
NLO corrections are so large that dyN{r)/N{r) is nega¬ 
tive around the “front” r ~ 1/Qs) which makes the so¬ 
lution progress unphysically in the wrong direction, with 
Qs deceasing with rapidity. For larger Qg, the behavior 
around r ~ l/Qs is less problematic, and we can focus 
on the small r tail of the amplitude. Here note that 
if dyN{r)/N(r) has a constant positive value, the ampli¬ 
tude grows exponentially in rapidity, but retains its shape 
as a function of r, resembling the small r behavior of the 
leading order evolution equation. This is indeed what 
happens for 7 = 0.6 and, marginally, for 7 = 0.8. For 
7 = 1 . 0 , however, we observe a negative, logarithmically 
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Q.s.o/Aqcd = 19 


<5.s,o/Aqcd = 19 


<5.v,o/Aqcd = 19 





Figure 7: Anomalous dimension 7 (r) = dln7V(r)/dlnr^ as a function of dipole size at different rapidities. The plots from right 
to left are for different rapidities y = 1,5,30. Shown in each plot are the solutions to the evolution equation with different 
initial anomalous dimensions. The solid lines are the solutions of the LO equation and the dashed lines the non-conformal 
dipoles, with different dashed lines corresponding to different initial anomalous dimensions. 


Q.s.o/Aqcd = 19 


<5.s,o/Aqcd = 19 


<5.v,o/Aqcd = 19 





Figure 8: Same as Fig. 7 for the conformal dipoles. 


decreasing dyN{r)/N{r) ~ Inr for r —>■ 0. This means 
that the evolution drives the amplitude towards a steeper 
shape, which in turn causes dyN(r)/N(r) to become even 
more negative for small r. Eventually this unstable de¬ 
velopment leads to a singularity in dyN{r)/N(r), which 
means that N{r) would cross zero at a finite r. At this 
point the integral on the r.h.s. of the BK equation would 
become divergent, so we impose N(r) > 0 by hand. The 
same evolution speed at y = 5 is shown in Fig. 4, where 
it can be seen that the evolution speed remains sensitive 
to the details of the initial condition even at large rapidi¬ 
ties. The corresponding results for the conformal dipole 
are shown in Figs. 5 and 6 , which show that the evolution 
speed of the conformal dipole is equally sensitive to the 
details of the initial condition. 

To understand better this unstable behaviour of the 
dipole amplitude shape we calculate the anomalous di¬ 
mension, defined as a logarithmic derivative of the dipole 
amplitude. 


7(r) 


din A^(r) 
dlnr^ 


( 12 ) 


The anomalous dimensions for the conformal and non- 


conformal dipoles at different rapidities are compared 
with the leading order BK solution in Figs. 7 and 8 . The 
instability of the solution is clearly visible at larger initial 
anomalous dimensions: with 7 = 1 in the initial condi¬ 
tion the anomalous dimension at small dipoles grows very 
large already after a few rapidity steps. With 7 = 0.8, 
a significantly longer evolution is needed before the so¬ 
lution becomes unstable. When the initial anomalous 
dimension is smaller (here 7 = 0 . 6 ), the unstable region 
is not reached within the rapidity interval studied here. 
Note that the solution with 7=1 does not evolve signif¬ 
icantly from 2 / = 5 to y = 30, as the evolution is domi¬ 
nated by distance scales where the the dipole amplitude 
would have already become negative in the numerical so¬ 
lution and has been frozen to iV(r) = 0. 

Contributions of different terms in the evolution equa¬ 
tion (1) are shown in Fig. 9 for the initial condition. In 
the figure the next-to-leading order contributions com¬ 
ing from the double logarithmic ^ Og In InK^/r^ 

and from the other ^ terms are shown separately. If 
the anomalous dimension in the initial condition is large 
enough, the double logarithmic term drives the evolution 
speed and is responsible for eventually turning the ampli- 
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Q.s.o/Aqcd ^ 19 Q,s,o/Aqcd ^ 19 




13,s,o/Aqcd ^ 19 



Figure 9: Evolution speed of the dipole amplitude at the initial condition. Shown are separately the full NLO and LO evolution 
equation results and the contributions from the conformal (no double logarithmic term) and non-conformal (only double 
logarithmic term) parts of the NLO BK equation. To demonstrate the location of the saturation scale the dipole amplitude is 
also shown as a thick line on the background. 





(a) 7 = 0.6 (b) 7 = 0.8 (c) 7 = 1.0 


Figure 10: Evolution speed of the conformal dipole amplitude at the initial condition. Shown are separately the full NLO and 
LO evolution equation results and the contributions from the In term and the rest of the evolution equation. To demonstrate 
the location of the saturation scale the dipole amplitude is also shown as a thick line on the background. 


tude negative at small dipoles. With small 7 both NLO 
contributions approach zero in the small dipole limit. For 
large dipoles in the saturation regime the NLO contribu¬ 
tions are still negative, but not large enough to change 
the sign of the evolution speed for Qs,o/Aqcd = 19. For 
a phenomenologically more relevant smaller value of Qs,o 
the coupling is larger and the NLO corrections make the 
evolution speed negative for all dipole sizes, as seen in 
Figs. 3 and 5. 

In the evolution equation of the conformal dipole, 
Eq. (7), the double logarithmic term is absent, but the 
definition of the conformal dipole makes an additional 
term ^ Oglnr^ appear in the NLO evolution equation. 
It is now exactly this In term which drives the con¬ 
formal dipole negative at small dipoles if the anoma¬ 
lous dimension at the initial condition is large enough. 
This is clearly seen in Fig. 10, where the NLO contribu¬ 
tions coming from the ^ In part and from the other 
NLO terms are shown. Note that in the original NLO 
BK equation the only logarithm in addition to the prob¬ 
lematic double logarithmic term is lnX^F'^/(X'^K^), 


which vanishes at r = 0. Compared to the evolution 
of the “non-conformal” dipole, the total evolution speed 
becomes negative at significantly smaller dipoles. If a 
smaller anomalous dimension is used at the initial condi¬ 
tion, also the contribution from the Inr^ term vanishes 
for small dipoles. 


In Ref. [32] an evolution equation for the conformal 
dipole in iV = 4 SYM theory is derived. We have 
checked that the conformal dipole in TV = 4 SYM has 
the same characteristic features as what was above shown 
in the case of QCD. Similarly we have also solved the 
evolution equation without using the Balitsky prescrip¬ 
tion for the running coupling in kernel Ki but instead 
absorbing terms proportional to /3 in the definition of 
as(Ai^) and using the smallest dipole prescription for as 
by choosing the scale to be set by the smallest dipole, 
1 ^^ The character¬ 

istic features of the solutions do not change in this change 
of the running coupling prescription. 










































V. CONCLUSIONS 

We have presented the first numerical solution to the 
next to leading order Balitsky-Kovchegov equation. The 
NLO corrections are shown to decrease the evolution 
speed and to be sensitive on the details of the initial 
condition. The slower evolution speed obtained by solv¬ 
ing the NLO evolution equations compared to the leading 
order BK is anticipated; in LO BK fits to HERA data 
the evolution speed needs to be reduced by evaluating 
the running coupling at a higher scale; see discussion in 
Refs. [3, 33]. However, as long as the solution to the evo¬ 
lution equation can become unphysical, too strong con¬ 
clusions on the effect of the higher order corrections on 
the evolution speed can not be made. 

The fact that the dipole amplitude may, depending on 
the initial condition, become negative and non vanishing 
at small dipoles is unphysical. This problem is only par¬ 
tially cured by writing the equation in terms of the con¬ 
formal dipole when the dipole amplitude becomes nega¬ 
tive only at significantly smaller dipoles. Even though it 
is possible to obtain an evolution which satisfies the re¬ 
quirement of a vanishing dipole amplitude at zero dipole 
size limit, we would like to get a stable evolution also in 


the case Qs,o ~ 1 GeV, 7 ~ 1 which has been the relevant 
region for phenomenological fits using the leading order 
equation (note that the leading order fits prefer values 
7 > 1 [1]). We have shown that the problematic behav¬ 
ior of the equation is associated with large logarithms of 
the small parent dipole size r, which corresponds to large 
transverse momentum. This confirms the result of [25] 
and suggests that the BK equation would indeed require 
a resummation of the same contributions that have been 
discussed in the context of the NLO BFKL equation [26- 
29]. This calls for a better understanding of the NLO 
evolution equation before the NLO dipole amplitude can 
be used in phenomenological applications. 
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